Metagenome Analysis

Only ancestral, streptomycin and control treatments are included, unless otherwise stated.

Authors

Niina Smolander

Shane Hogle

Libraries and functions

Prior to this analysis, the variants have been filtered using Mutect2’s FilterMutectCalls as well as separately for SNPs and INDELs using the following hard filtering factors: SOR, FS, MQ, MQRankSum, ReadPosRankSum, SNP/INDEL clustering (3/35) and SNP/INDEL proximity.

Code
library(tidyverse)
library(Polychrome)
library(withr)
library(fs)
library(compositions)


plot_nuc_survival <- function(df, ylim, ybreaks, ncol){
  ggplot(df, aes(x=m, y=value, color=name, linetype = name)) + 
    geom_step(direction = "mid") +
    scale_x_continuous(limits = c(0.5, 3.5), 
                       breaks = c(1, 2, 3), 
                       minor_breaks = NULL) +
    scale_y_log10(breaks = {{ ybreaks }},
                  labels = scales::trans_format("log10", scales::math_format(10^.x))) + 
    labs(x=("Nucleotide multiplicitiy, m"), 
         y=("Total mutations ≥ m"),
         color = NULL) +
    scale_color_manual(values = c("blue", "red")) +
    scale_linetype(guide = 'none') +
    facet_wrap(~label, ncol = {{ ncol }}) +
    annotation_logticks(sides = "l", color = "grey40") +
    coord_cartesian(ylim = c(1, {{ ylim }})) +
    theme_bw() +
    theme(
      strip.placement = 'outside',
      strip.background = element_blank(),
      panel.grid = element_blank(),
      legend.position = "bottom"
    )
}

Variants and annotations

Code
varfiles <- list.files(path = "../Amanda_metagenomedata/tables_better_annotated/")

varfiles_snpeff <- varfiles[base::grep("snpeff", varfiles)]

variants <- list()
for(i in varfiles_snpeff){
  name <- str_remove(i, "_L003.VF.HF.snpeff.tsv")
  variants[[name]] <- read_tsv(paste0("../Amanda_metagenomedata/tables_better_annotated/", i),
                               skip = 1,
                               col_names = c("CHROM",
                                             "POS",
                                             "REF",
                                             "ALT",
                                             "FILTER",
                                             "EFFECT",
                                             "IMPACT",
                                             "GENE",
                                             "BIOTYPE",
                                             "HGVS_C",
                                             "HGVS_P",
                                             "CDS_POS",
                                             "CDS_LEN",
                                             "AA_POS",
                                             "AA_LEN"),
                               col_types = c("cdcccccccccdddd"))
}


variants <- bind_rows(variants, .id = "Sample") %>%
  filter(FILTER == "PASS")

varfiles_variants <- varfiles[base::grep("variants", varfiles)]

variants2 <- list()
for(i in varfiles_variants){
  name <- str_remove(i, "_L003.VF.HF.variants.tsv")
  variants2[[name]] <- read_tsv(paste0("../Amanda_metagenomedata/tables_better_annotated/", i),
                                skip = 1,
                                col_names = c("CHROM",
                                              "POS",
                                              "REF",
                                              "ALT",
                                              "FILTER",
                                              "TYPE",
                                              "AD",
                                              "AF",
                                              "DP",
                                              "FAD",
                                              "GQ",
                                              "VAF"),
                                col_types = c("cdcccccccccc"))
}


variants2 <- bind_rows(variants2, .id = "Sample") %>%
  filter(FILTER == "PASS")

# Combine and format rows with multiple ALTs

variants_all <- full_join(variants, variants2,
                          by = join_by(Sample, CHROM, POS, REF, ALT, FILTER)) %>%
  separate_wider_delim(cols = "Sample", delim = "-", names = c("HAMBI", "Sample"))


variants_all_tidy <- variants_all %>%
  separate_wider_delim(cols = "AD", delim = ",", names = c("REF_AD" ,"ALT1_AD", "ALT2_AD", "ALT3_AD", "ALT4_AD"), too_few = "align_start") %>%
  separate_wider_delim(cols = "FAD", delim = ",", names = c("REF_FAD" ,"ALT1_FAD", "ALT2_FAD", "ALT3_FAD", "ALT4_FAD"), too_few = "align_start") %>%
  mutate(ALT_AD = paste(ALT1_AD, ALT2_AD, ALT3_AD, ALT4_AD, sep = ","),
         ALT_FAD = paste(ALT1_FAD, ALT2_FAD, ALT3_FAD, ALT4_FAD, sep = ",")) %>%
  mutate(ALT_AD = str_remove_all(ALT_AD, ",NA"),
         ALT_FAD = str_remove_all(ALT_FAD, ",NA")) %>%
  select(-c("ALT1_AD", "ALT2_AD", "ALT3_AD", "ALT4_AD", "ALT1_FAD", "ALT2_FAD", "ALT3_FAD", "ALT4_FAD", "FILTER", "GQ")) %>%
  separate_longer_delim(cols = c("ALT", "AF", "VAF", "ALT_AD", "ALT_FAD"), delim = ",") %>%
  mutate(across(c("AF", "DP", "VAF", "REF_AD", "ALT_AD", "REF_FAD", "ALT_FAD"), as.numeric))

Filtering

Keeping only the variants where allelic depth (AD) ≥ 5, coverage (DP) ≥ 5 and VAF ≥ 0.05. Also variants that are in mobile genetic elements are removed.

Code
# AD, DP and VAF
variants_all_tidy <- variants_all_tidy %>%
  filter(ALT_AD >= 5 & DP >= 5 & VAF >= 0.05)

# MGE
genomad <- read_tsv("../Amanda_metagenomedata/MGE/genomad_HAMBI_combined.tsv",
                    skip = 1,
                    col_names = c("HAMBI",
                                  "CHROM",
                                  "Start",
                                  "End",
                                  "Name",
                                  "prediction_tool"),
                    col_types = c("ccddcc"))

mgefinder <- read_tsv("../Amanda_metagenomedata/MGE/MGE_finder_HAMBI_combined.tsv",
                      skip = 1,
                      col_names = c("HAMBI",
                                    "Name",
                                    "CHROM",
                                    "Start",
                                    "End",
                                    "prediction_tool"),
                      col_types = c("cccddc"))


variants_all_tidy_mge <- variants_all_tidy %>% 
  left_join(mgefinder, by = join_by(HAMBI, CHROM), relationship = "many-to-many") %>% 
  filter(!(POS >= Start & POS <= End)) %>% 
  dplyr::select(-Name, -Start, -End, -prediction_tool) %>% 
  distinct() %>% 
  left_join(genomad, by = join_by(HAMBI, CHROM), relationship = "many-to-many") %>% 
  filter(!(POS >= Start & POS <= End)) %>% 
  dplyr::select(-Name, -Start, -End, -prediction_tool) %>% 
  distinct()

# metadata
meta <- read_csv("../Amanda_metagenomedata/Experiment_info/key.csv") %>%
  mutate(sample = str_remove_all(sample, "_L003"))

variants_all_tidy_final <- left_join(variants_all_tidy_mge, meta,
                                     by = join_by(Sample == sample)) %>%
  mutate(Experiment = ifelse(is.na(Experiment),
                             "t0",
                             Experiment))

# remove tetrahymena
variants_strep <- variants_all_tidy_final %>%
  filter(Experiment != "Tetrahymena" & Experiment != "Both") %>%
  filter(History != "Tetrahymena" & History != "Both")

Distance matrices using variants

Using Bray-Curtis dissimilarity.

Code
# prep data
var_matrix_all <- variants_all_tidy_final %>%
  mutate(var_id = paste(CHROM, POS, ALT, sep = "_")) %>%
  select(var_id, Experiment, History, Timepoint, Replicate, Sample, VAF) %>%
  pivot_wider(names_from = var_id, values_from = VAF, values_fill = 0) %>%
  as.data.frame(.)

rownames(var_matrix_all) <- var_matrix_all$Sample

var_matrix <- variants_strep %>%
  mutate(var_id = paste(CHROM, POS, ALT, sep = "_")) %>%
  select(var_id, Experiment, History, Timepoint, Replicate, Sample, VAF) %>%
  pivot_wider(names_from = var_id, values_from = VAF, values_fill = 0) %>%
  as.data.frame(.)

rownames(var_matrix) <- var_matrix$Sample

# capscale
cap_results_all <- vegan::capscale(var_matrix_all[-(1:5)] ~ History*Experiment,
                               data = var_matrix_all,
                               dist = "bray")


cap_results <- vegan::capscale(var_matrix[-(1:5)] ~ History*Experiment,
                               data = var_matrix,
                               dist = "bray")

# plot prep all
baseplot_all <- plot(cap_results_all)
Code
sites_all <- data.frame(baseplot_all$sites)
sites_all$Sample <- rownames(sites_all)
sites_all <- dplyr::left_join(sites_all, select(var_matrix_all, History, Experiment, Timepoint, Replicate, Sample), by = join_by(Sample))

# plot prep strep
baseplot <- plot(cap_results)
Code
sites <- data.frame(baseplot$sites)
sites$Sample <- rownames(sites)
sites <- dplyr::left_join(sites, select(var_matrix, History, Experiment, Timepoint, Replicate, Sample), by = join_by(Sample))

CAP plot for the whole data (incl. tetrahymena samples)

First plot split by History, the second one split by Experiment. Experiment t0 refers to the timepoint 0 (start of the experiment).

Code
ggplot() +
  geom_point(data=sites_all, aes(x=CAP1, y=CAP2, color=Experiment), alpha = 1.5, size = 4) +
  labs(x = paste0("CAP1 (", round(summary(vegan::eigenvals(cap_results_all))[2,1] *100, 2)
, "%)"),
       y = paste0("CAP2 (", round(summary(vegan::eigenvals(cap_results_all))[2,2] *100, 2)
, "%)")) +
  facet_wrap(History ~ .)

Code
ggplot() +
  geom_point(data=sites_all, aes(x=CAP1, y=CAP2, color=History), alpha = 1.5, size = 4) +
  labs(x = paste0("CAP1 (", round(summary(vegan::eigenvals(cap_results_all))[2,1] *100, 2)
, "%)"),
       y = paste0("CAP2 (", round(summary(vegan::eigenvals(cap_results_all))[2,2] *100, 2)
, "%)")) +
  facet_wrap(Experiment ~ .)

CAP plots for the no-tetrahymena data

First plot split by History, second one split by Experiment. Experiment t0 refers to the timepoint 0 (start of the experiment).

Code
ggplot() +
  geom_point(data=sites, aes(x=CAP1, y=CAP2, color=Experiment), alpha = 1.5, size = 4) +
  labs(x = paste0("CAP1 (", round(summary(vegan::eigenvals(cap_results))[2,1] *100, 2)
, "%)"),
       y = paste0("CAP2 (", round(summary(vegan::eigenvals(cap_results))[2,2] *100, 2)
, "%)")) +
  facet_wrap(History ~ .)

Code
ggplot() +
  geom_point(data=sites, aes(x=CAP1, y=CAP2, color=History), alpha = 1.5, size = 4) +
  labs(x = paste0("CAP1 (", round(summary(vegan::eigenvals(cap_results))[2,1] *100, 2)
, "%)"),
       y = paste0("CAP2 (", round(summary(vegan::eigenvals(cap_results))[2,2] *100, 2)
, "%)")) +
  facet_wrap(Experiment ~ .)

Distances between the timepoints t3 and t0

Distance between the t3 and t0 points in the CAP plot for each History-Experiment combination where two timepoints existed.

Code
test_t0 <- sites %>% select(-Sample) %>%
  filter(Timepoint == "t0")

test_t3 <- sites %>% select(-Sample) %>%
  filter(Timepoint == "t3" & History != "Ancestor")

test <- right_join(test_t0, test_t3,
          by = join_by(History, Replicate),
          suffix = c("_t0", "_t3")) %>%
  select(-c("Experiment_t0", "Timepoint_t0", "Timepoint_t3")) %>%
  rename("Experiment" = "Experiment_t3") %>%
  mutate(CAP1_dist = (CAP1_t0 - CAP1_t3)^2,
         CAP2_dist = (CAP2_t0 - CAP2_t3)^2) %>%
  mutate(dist = sqrt(CAP1_dist + CAP2_dist)) %>%
  mutate(History_Experiment = paste(History, Experiment, sep = "_"))

test %>%
  ggplot(aes(x = History_Experiment, y = dist)) +
  geom_boxplot() +
  labs(y = "Distance between t3 and t0")

Parallelism at nucleotide level

Looking into how many mutations are found in how many replicates in different HAMBIs at the timepoint 3. Fixed variants (VAF = 1) seen in the timepoint 0 are excluded (but not replicate by replicate). Using https://slhogle.github.io/hambiDoubleEvo/R/metagenome/02_parallelism.html as reference.

Code
chrom_lengths <- read_tsv("../Amanda_metagenomedata/Experiment_info/chrom_lengths.tsv")

# exclude fixed variants from timepoint 0
tofilter <- variants_strep %>% 
  filter(Timepoint == "t0") %>% 
  filter(VAF == 1) %>% 
  dplyr::select(CHROM, POS, REF, ALT, History) %>% 
  distinct()


variants_filt <- anti_join(variants_strep, tofilter, by = join_by(CHROM, POS, REF, ALT, History)) %>%
  filter(Timepoint != "t0")

HAMBI_mutations_filt <- variants_filt %>%
  group_by(HAMBI, Experiment, History, CHROM, POS, REF, ALT) %>% 
  mutate(mi=n_distinct(Replicate)) %>% ungroup() %>%
  group_by(HAMBI) %>% 
  count(mi) %>% ungroup() %>%
  group_by(mi) %>% 
  mutate(tot = sum(n)) %>% 
  arrange(desc(mi)) %>% ungroup()

HAMBI_mutations_filt2 <- variants_filt %>%
  group_by(HAMBI, Experiment, History, CHROM, POS, REF, ALT) %>% 
  mutate(mi=n_distinct(Replicate)) %>% ungroup() %>%
  group_by(HAMBI, History, Experiment) %>% 
  count(mi) %>% ungroup() %>%
  group_by(mi) %>% 
  mutate(tot = sum(n)) %>% 
  arrange(desc(mi)) %>% ungroup()

Table

mi: number of replicates.

Code
knitr::kable(HAMBI_mutations_filt)
HAMBI mi n tot
HAMBI_0006 3 63 114
HAMBI_1842 3 9 114
HAMBI_1972 3 30 114
HAMBI_1977 3 3 114
HAMBI_2659 3 9 114
HAMBI_0006 2 22 92
HAMBI_0403 2 2 92
HAMBI_1842 2 8 92
HAMBI_1972 2 10 92
HAMBI_1977 2 28 92
HAMBI_2659 2 22 92
HAMBI_0006 1 36 309
HAMBI_0403 1 10 309
HAMBI_1287 1 5 309
HAMBI_1842 1 38 309
HAMBI_1972 1 31 309
HAMBI_1977 1 77 309
HAMBI_2164 1 1 309
HAMBI_2443 1 11 309
HAMBI_2659 1 100 309

Plot

Distribution of nucleotide multiplicity (the number of replicates with the same genomic mutation) for HAMBI x history x treatment combinations. The observed data in blue, the null expectation in red.

Code
nuc_survival <- variants_filt %>%
  dplyr::select(HAMBI, History, Experiment, POS, REF, ALT, Replicate) %>% 
  # multiplicity = number of replicate populations each unique mutation is
  # observed in
  summarize(m = n_distinct(Replicate), .by = c(HAMBI, History, Experiment, POS, REF, ALT)) %>%
  # now calculate the total number of mutations across all replicates so we need
  # to ungroup by mutation position/alt allele but because we still want to
  # determine this value by treatment category we keep the treatment category
  # grouping. However, this should be changed if you want to for example average
  # over all the treatment conditions on a species basis
  group_by(HAMBI, History, Experiment) %>%
  count(m, name = "m_count") %>%
  mutate(n = m * m_count,
         Ntot = sum(n),
         perc = n / Ntot * 100) %>%
  left_join(chrom_lengths, by = join_by(HAMBI)) %>%
  arrange(cur_group_id(), desc(m)) %>%
  #  dpois() tells the probability mass at a given number of counts. Here we
  #  want to get the probability of observing n mutations with multiplicity
  #  = mi (i.e. the counts of mi in the observed data). We assume that
  #  mutations independently occur on the genome of size Ltot at a rate of
  #  lambda = Ntot/Ltot and that generally the events are rare. Thus this
  #  situation can be modeled by the Poisson distribution. We can get the
  #  binned number of mutations per level of multiplicity m by multiplying
  #  the probability by the length of the genome and the binned mutations
  #  divided by the total number of mutations.
  mutate(m_count_expected = cumsum((m_count / Ntot) *
                                     sum_bases *
                                     dpois(m, lambda = Ntot / sum_bases))) %>%
  #dplyr::select(-num_contigs) %>%
  relocate(m, n, Ntot, perc, m_count, m_count_expected) %>%
  ungroup()

# setup for plotting
nuc_survival_plot <- nuc_survival %>% 
  group_by(HAMBI, History, Experiment) %>% 
  # when there is only one multiplicity observed for a mutation filter such
  # that the multiplicty of that mutation must be greater than one.
  # Otherwise include all remaining mutations (m > 0)
  filter(case_when(n() == 1 ~ m > 1,
                   TRUE ~ m > 0)) %>% 
  pivot_longer(cols = c("m_count", "m_count_expected")) %>% 
  mutate(label = paste(HAMBI, "History:", History, "\nExperiment:", Experiment)) %>% 
  # and make final plot
  plot_nuc_survival(., 5000, c(1, 10, 100, 1000, 10000), 4)

nuc_survival_plot

Barplot of the mutation distributions

Split by history and experiment.

Code
HAMBI_mutations_filt2 %>%
  ggplot(aes(x = mi, y = n, fill = HAMBI)) +
  geom_bar(stat = "identity") +
  facet_grid(Experiment ~ History)

Table

Distributions summarised.

Code
nuc_survival_table <- nuc_survival %>% 
  mutate(Ntot_f = sum(n)) %>%
  reframe(perc = sum(n)/Ntot_f*100, .by = "m") %>% 
  distinct() %>% 
  mutate(cum_perc = cumsum(perc))

knitr::kable(nuc_survival_table)
m perc cum_perc
3 22.13592 22.13592
2 17.86408 40.00000
1 60.00000 100.00000

Variants that existed in YSK

Plot of variants that existed in t0 (YSK) samples in at least two replicates. Plot split by experimental condition, NA panel has the variants that were seen in YSK but not in timepoint 3 (t3). Color indicates the difference in variant allele frequency between t3 and t0 (t3 minus t0). Fixed variants (VAF >= 0.9) marked with the text “YSK” / “t3” / “Both”.

Code
gene_name_key <- read_tsv("../Amanda_metagenomedata/Experiment_info/gene_names_and_products_key.tsv") %>%
  select(ID, gene_name)

variants_strep <- variants_strep %>%
  mutate(GENE1 = ifelse(str_detect(GENE, "-"),
                        str_split_i(GENE, "-", 1),
                        NA),
         GENE2 = ifelse(str_detect(GENE, "-"),
                        str_split_i(GENE, "-", 2),
                        NA)) %>%
  left_join(gene_name_key, by = join_by(GENE1 == ID)) %>%
  rename(gene_name1 = gene_name) %>%
  left_join(gene_name_key, by = join_by(GENE2 == ID)) %>%
  rename(gene_name2 = gene_name) %>%
  mutate(gene_name1 = ifelse(is.na(gene_name1),
                             "?",
                             gene_name1),
         gene_name2 = ifelse(is.na(gene_name2),
                             "?",
                             gene_name2)) %>%
  left_join(gene_name_key, by = join_by(GENE == ID)) %>%
  mutate(gene_name = ifelse(str_detect(GENE, "-"),
                            paste(gene_name1, gene_name2, sep = " - "),
                            gene_name),
         gene_name = ifelse(is.na(gene_name),
                            "?",
                            gene_name))

# collect which replicates in t3
t3_samples <- variants_strep %>% filter(Timepoint == "t3") %>%
  mutate(id = paste(History, Replicate, sep = "_")) %>%
  pull(id) %>% unique()

t3_variants <- variants_strep %>% filter(Timepoint == "t3") %>%
  select(HAMBI, CHROM, POS, REF, ALT, GENE, EFFECT, History, Experiment, Replicate, VAF, gene_name) %>%
  rename("VAF_t3" = VAF)

YSK_variants <- variants_strep %>%
  filter(Timepoint == "t0") %>%
  select(HAMBI, CHROM, POS, REF, ALT, GENE, gene_name, EFFECT, History, Replicate, VAF) %>%
  rename("VAF_t0" = VAF) %>%
  filter(paste(History, Replicate, sep = "_") %in% t3_samples) %>%
  group_by(CHROM, POS, ALT, History) %>%
  mutate(mi = n_distinct(Replicate)) %>% ungroup() %>%
  left_join(t3_variants,
            by = join_by(HAMBI, CHROM, POS, REF, ALT, GENE, EFFECT, History, Replicate, gene_name)) %>%
  mutate(VAF_t3 = ifelse(is.na(VAF_t3), 0, VAF_t3),
         VAF_diff = VAF_t3 - VAF_t0)

YSK_variants %>% filter(mi >= 2) %>%
  mutate(ALT = ifelse(ALT == "CGCCGGCAAGCCAGCTCCCGCAGGTACTGCACAAGCCTTGAAGGCAGTGATGTCCCTGTGGGAGCTGGCTTGCCGGCGATAGG",
                      "long_ins_1",
                      ALT),
         GENE_ALT = paste0(GENE,"_", POS, "_", ALT, " ( ", gene_name, " ) "),
         History_Replicate = paste(str_sub(History, 1, 5), Replicate, sep = "_"),
         Fixed = ifelse(VAF_t0 >= 0.9 & VAF_t3 >= 0.9,
                        "Both",
                        ifelse(VAF_t0 >= 0.9 & VAF_t3 < 0.9,
                               "YSK",
                               ifelse(VAF_t0 < 0.9 & VAF_t3 >= 0.9,
                                      "t3",
                                      "")))) %>%
  ggplot(aes(x = History_Replicate, y = GENE_ALT, fill = VAF_diff)) +
  geom_tile() +
  geom_text(aes(label = Fixed), color = "white") +
  #scale_fill_viridis_c(option = "D", direction = 1) +
  scale_fill_gradient2(low = "darkblue",
                       mid = "grey",
                       high = "green") +
  facet_grid(Experiment ~ ., scales = "free", space = "free")

De novo variants found in all three replicates

Code
t3_variants_mi3 <- t3_variants %>%
  anti_join(YSK_variants,
            by = join_by(HAMBI, CHROM, POS, REF, ALT, GENE, History, Replicate)) %>%
  group_by(CHROM, POS, ALT, History, Experiment) %>%
  mutate(mi = n_distinct(Replicate)) %>% ungroup() %>%
  filter(mi >= 3)
  

t3_variants_mi3 %>%
  mutate(GENE_ALT = paste0(GENE,"_", POS, "_", ALT, " ( ", gene_name, " ) "),
         History_Replicate = paste(str_sub(History, 1, 5), Replicate, sep = "_"),
         Fixed = ifelse(VAF_t3 >= 0.9,
                        "Fixed",
                        "")) %>%
  ggplot(aes(x = History_Replicate, y = GENE_ALT, fill = VAF_t3)) +
  geom_tile() +
  geom_text(aes(label = Fixed), color = "white") +
  scale_fill_viridis_c(option = "D", direction = 1) +
  # scale_fill_gradient2(low = "darkblue",
  #                      mid = "grey",
  #                      high = "green") +
  facet_grid(Experiment ~ ., scales = "free", space = "free")

De novo variants found in two replicates

Code
t3_variants_mi2 <- t3_variants %>%
  anti_join(YSK_variants,
            by = join_by(HAMBI, CHROM, POS, REF, ALT, GENE, History, Replicate)) %>%
  group_by(CHROM, POS, ALT, History, Experiment) %>%
  mutate(mi = n_distinct(Replicate)) %>% ungroup() %>%
  filter(mi == 2)

t3_variants_mi2 %>%
  mutate(ALT = ifelse(ALT == "CGCCGGCAAGCCAGCTCCCGCAGGTACTGCACAAGCCTTGAAGGCAGTGATGTCCCTGTGGGAGCTGGCTTGCCGGCGATAGG",
                      "long_ins_1",
                      ALT),
         GENE_ALT = paste0(GENE,"_", POS, "_", ALT, " ( ", gene_name, " ) "),
         History_Replicate = paste(str_sub(History, 1, 5), Replicate, sep = "_"),
         Fixed = ifelse(VAF_t3 >= 0.9,
                        "Fixed",
                        "")) %>%
  ggplot(aes(x = History_Replicate, y = GENE_ALT, fill = VAF_t3)) +
  geom_tile() +
  geom_text(aes(label = Fixed), color = "white") +
  scale_fill_viridis_c(option = "D", direction = 1) +
  # scale_fill_gradient2(low = "darkblue",
  #                      mid = "grey",
  #                      high = "green") +
  facet_grid(Experiment ~ ., scales = "free", space = "free")